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Abstract 

Images  are  used  for  both  aerial  and  space  imagery  applications,  including  target 
detection  and  tracking.  The  current  problem  concerning  objects  in  geosynchronous 
orbit  is  that  they  are  dim  and  hard  to  resolve  because  of  their  distance.  This  work  will 
further  the  combined  effort  of  AFIT  and  AFRL  to  provide  enhanced  space  situational 
awareness  (SSA)  and  space  surveillance.  SSA  is  critical  in  a  time  when  many  countries 
possess  the  technology  to  put  satellites  into  orbit.  Enhanced  imaging  technology 
improves  the  Air  Force’s  ability  to  see  if  foreign  satellites  or  other  space  hardware  are 
operating  in  the  vicinity  of  our  own  assets  at  geosynchronous  orbit.  Image  deblurring 
or  denoising  is  a  crucial  part  of  restoring  images  that  have  been  distorted  either 
by  movement  during  the  capture  process,  using  out-of-focus  optics,  or  atmospheric 
turbulence.  The  goal  of  this  work  is  to  develop  a  new  blind  deconvolution  method  for 
imaging  objects  at  geosynchronous  orbit.  It  will  feature  an  expectation  maximization 
(EM)  approach  to  iteratively  deblur  an  image  while  using  the  convergence  of  the 
image’s  variance  as  the  stopping  criteria 
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Blind  Deconvolution  Method  Oe  Image  Deblurring  Using 

Convergence  Oe  Variance 

I.  Introduction 

This  chapter  describes  the  problem  to  be  addressed  by  this  research.  Background  on 
the  problem  and  research  goals  are  also  provided. 

1 . 1  Motivation 

Space  Situational  Awareness  (SSA)  is  a  key  mission  of  the  United  States  Air 
Force  Space  Command.  One  aspect  of  SSA  involves  using  both  telescope  networks 
and  radars  to  detect,  identify,  record  and  track  all  man-made  objects  orbiting  the 
earth.  Knowing  the  exact  locations  of  these  orbiting  objects  in  space  is  crucial  for 
future  space  operation  safety.  The  SSA  mission  has  become  even  more  important 
with  recent  events  such  as  Iridium  and  Cosmos  satellite  collision  and  the  Chinese 
anti-satellite  missile  test  in  2007,  both  of  which  created  large  swaths  of  space  debris. 
This  debris,  known  as  space  junk,  will  be  an  ongoing  risk  to  US  satellites  for  years  to 
come  as  the  orbits  of  the  debris  degrade  toward  Earth. 

1.2  Background 

Images  are  used  for  both  aerial  and  space  imagery  applications,  including  target 
detection  and  tracking.  The  current  problem  concerning  objects  in  geosynchronous 
orbit  is  that  they  are  dim  and  hard  to  resolve  because  of  their  distance.  This  work  will 
further  the  combined  effort  of  AFIT  and  AFRL  to  provide  enhanced  SSA  and  space 
surveillance.  SSA  is  critical  in  a  time  when  many  conntries  possess  the  technology  to 
put  satellites  into  orbit  and  enhanced  imaging  technology  will  improve  the  Air  Force’s 
ability  to  see  if  foreign  satellites  or  other  space  hardware  are  operating  in  the  vicinity 
of  our  own  assets  at  geosynchronous  orbit. 


1 


Image  deblurring  or  denoising  is  a  crucial  part  of  restoring  images  that  have  been 
distorted  either  by  movement  during  the  capture  process,  using  out-of-focus  optics, 
or  atmospheric  turbulence.  There  currently  exist  a  number  of  methods  for  image 
deblurring.  These  either  use  the  short  optical  transfer  function  (OTF),  which  negates 
atmospheric  tilt  caused  by  phase  delay,  or  rely  on  a  posteriori  data  in  order  to  make 
an  estimate  abont  the  seeing  parameter,  which  describes  the  effect  of  atmospheric 
tnrbnlence.  The  goal  of  this  work  is  to  develop  a  new  blind  deconvolntion  method  for 
long  exposnre  imaging  of  objects  at  geosynchronons  orbit  withont  the  need  for  any  a 
posteriori  data. 

The  proposed  research  effort  wonld  attempt  to  produce  a  blind  deconvolntion 
algorithm  that  prodnces  usefnl  results  when  processing  astronomical  scenes.  This  new 
techniqne  will  capitalize  on  the  statistics  of  the  blnrry  image  and  the  rehned  image 
estimate,  in  an  iterative  approach  to  converge  on  the  correct  seeing  parameter.  This 
approach  has  the  potential  to  perform  blind  deconvolntion  on  myriad  data  sonrces 
and  applications. 

1.3  Research  Goals 

The  expected  result  of  this  work  is  the  development  of  a  new  approach  for 
blind  deconvolntion  that  improves  on  cnrrent  methods  for  deblnrring  astronomical 
data.  Experiments  condncted  using  simulated  data  will  feature  binary  star  images  as 
well  as  images  of  satellites  blurred  by  atmospheric  turbulence.  The  signal  to  noise 
ratio  of  the  data  will  also  be  varied  via  changes  in  the  coherence  parameter  of  the 
simnlated  laser  illnmination.  Experiments  using  measnred  data  will  use  the  value  of 
the  estimated  seeing  parameter  as  a  metric  for  snccess  as  gronnd  trnth  measnrements 
of  that  parameter  are  available  for  each  set  of  measnred  imagery. 

1.4  Thesis  Organization 

Chapter  II  provides  a  description  of  the  Maximnm  A  Posteriori  (MAP)  es¬ 
timation  algorithm  and  related  work,  attempting  to  solve  the  blind  deconvolntion 
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problem  using  the  long  exposure  atmospheric  transfer  function.  Chapter  III  details 
the  methodology  used  in  this  research,  and  the  simulations  and  experiments  which 
were  conducted.  Chapter  IV  yields  the  results  from  the  simulations  and  experiments 
described  in  Chapter  III.  Finally,  Chapter  V  provides  a  summary  for  this  research 
and  suggests  potential  areas  for  follow-on  research. 
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II.  Background 

This  chapter  provides  the  technical  background  necessary  for  understanding  the  over¬ 
all  concepts  of  this  research.  First,  a  description  of  blind  deconvolution  is  provided. 
Next,  related  work  is  presented  including  an  examination  of  both  the  MAP  estimation 
and  APEX  techniques.  The  APEX  method,  which  is  not  an  acronym,  was  developed 
by  Alfred  Carasso  [5].  The  MAP  algorithm  is  provided  along  with  its  derivation. 
Examples  are  provided  to  demonstrate  how  these  algorithms  perform  on  fully  illu¬ 
minated  datasets.  Finally,  examples  are  given  to  show  how  these  algorithms  are  not 
useful  when  processing  atmospheric  data,  and  an  alternative  method  is  suggested. 

2.1  Blind  Deconvolution 

Deconvolution  is  the  process  of  hltering  a  signal  to  compensate  for  an  undesired 
convolution.  This  convolution  operation  on  the  data  is  caused  by  the  effect  of  atmo¬ 
spheric  turbulence.  The  goal  of  deconvolution  is  to  recreate  the  signal  as  it  existed 
before  the  convolution  took  place.  This  usually  requires  that  the  characteristics  of 
the  convolution  are  known.  In  this  case  it  is  assumed  that  the  point  spread  function 
(PSF)  responsible  for  the  blurring  is  known. 

Blind  deconvolution,  on  the  other  hand,  is  a  deconvolution  technique  that  per¬ 
mits  the  recovery  of  the  target  scene  from  either  a  single  or  a  set  of  blurred  images. 
This  can  still  be  achieved  with  a  poorly  determined  or  unknown  PSF.  In  blind  decon¬ 
volution  techniques,  the  PSF  is  estimated  from  the  image  or  image  set,  allowing  the 
deconvolution  to  be  performed. 

Blind  deconvolution  can  be  performed  iteratively,  where  each  iteration  improves 
the  estimation  of  the  PSF  and  the  scene,  as  in  the  case  of  MAP  Estimation.  It  can  also 
be  performed  non-iteratively  like  the  APEX  algorithm,  where  the  algorithm,  based  on 
exterior  information,  extracts  the  PSF.  Beginning  with  a  good  estimate  of  the  PSF 
is  helpful  for  quicker  convergence  but  is  not  necessary  to  recover  the  target  scene. 


4 


2.2  Optical  Transfer  Function 

2.2.1  Average  Short-Exposure  OTF.  The  average  short-exposure  OTF  can 
be  achieved  by  either  of  two  methods.  The  first  method  involves  averaging  several 
short-exposure  images  together.  The  second  requires  that  the  tilt  from  several  long- 
exposure  images  is  removed  before  they  are  averaged.  Short  exposure  images  are 
created  by  integrating  light  from  the  target  where  the  integration  time  does  not  exceed 
second  [8].  The  form  of  the  equation  for  the  short  exposure  transfer  function  can 
be  seen  here  in  (1). 


Hs{v)  =  exp 


-3.44 


Ro  ) 


5/3 


1  —  a 


(1) 


where  a  takes  on  the  value  unity  for  “near-held”  and  ^  for  “far-held”  propagation. 
A  is  the  wavelength  of  the  light,  v  is  its  frequency  and  /  is  the  focal  length  of  the 
imaging  system.  Rq  is  the  atmospheric  seeing  parameter  and  Dq  is  the  diameter  of 
the  optical  system.  When  a  is  set  to  zero,  the  equation  takes  on  the  form  of  the 
long-exposure  OTF,  in  (2). 

The  atmospheric  coherence  diameter  Rq,  also  known  as  the  seeing  parameter  hrst 
introduced  by  Fried  [8],  describes  the  image-degrading  ehects  of  the  atmosphere. 
Typical  values  of  Rq  at  a  mountain-top  observatory  might  range  from  5  cm  under 
poor  visibility  up  to  20  cm  under  exceptional  viewing  conditions. 


2.2.2  Long-Exposure  OTE.  Long  exposure  dihers  from  short  exposure  with 
respect  to  the  integration  time  used.  Long  exposure  images  typically  use  integration 
times  much  greater  than  second.  This  increased  integration  time  also  introduces 
wavefront  tilt  or  wavefront  fluctuations  caused  by  the  atmospheric  winds  [8].  The 
equation  for  the  long  exposure  transfer  function  can  be  seen  here  in  (2). 


Hl{v) 


(2) 
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2.3  Statistics 


The  two  statistics  that  are  necessary  for  this  thesis  are  the  mean  and  variance 
of  a  dataset.  The  mean  or  average  value  is  used  to  determine  the  signal  to  noise  ratio 
to  see  if  it  is  sufficient  for  this  method  to  be  able  to  converge  properly. 

2.3.1  Mean.  The  mean  of  a  random  variable  is  equal  to  its  average  or 
expected  value.  The  mean  of  the  images  used  in  this  thesis  is  derived  by  summing 
each  pixel  value  in  the  image  and  dividing  that  total  by  the  number  of  pixels  present. 


Mean{X)  =  E[X] 


i=l 


(3) 


where  X  represents  any  random  variable  and  E[-]  represents  the  expected  value  op¬ 
erator. 

2.3.2  Variance.  The  variance  of  a  random  variable  is  equal  to  the  mean 
of  the  squared  deviation  of  that  variable  from  its  expected  value.  The  variance  is, 
therefore,  a  measurement  of  how  much  fluctuation  exists  within  the  values  of  that 
variable.  Likewise,  the  variance  of  the  images  can  be  determined  by  hrst  hnding  the 
mean  or  average  pixel  value,  as  described  above.  Next  subtract  the  mean  from  each 
of  the  pixel  values  and  square  the  result.  Next  sum  up  all  these  results  and  divide  by 
the  number  of  pixels  to  determine  the  average  variance  of  the  scene. 
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N 


-  i-'f 


Var{X)  =  E[{X-E[x]f 


(4) 


where  X  represents  any  random  variable,  {X  —  E[x])  is  the  deviation  and  fi  is 
the  mean  as  in  (3). 

2.4  Signal  To  Noise  Ratio 

Signal  to  noise  ratio  (SNR)  refers  to  the  measure  of  a  signal’s  strength  relative 
to  any  background  noise.  If  the  SNR  in  question  is  not  particularly  strong  (i.e.  - 
orders  of  magnitude),  then  it  may  be  to  difficult  to  properly  isolate  the  signal  from 
the  noise  floor.  SNR  is  typically  represented  as  a  ratio  of  either  power  or  amplitude 
of  a  signal  with  respect  to  the  noise. 


SNR 


signal 


A 


signal 


jA-r. 


(5) 


where  Psignai  represents  the  power  of  the  signal,  Pnoise  is  the  power  of  the  noise  level, 
^signal  IS  the  amplitude  of  the  signal,  and  A^oise  is  the  amplitude  of  the  noise  level. 
Throughout  this  research  and  the  simulations,  SNR  will  be  referred  to  as  a  ratio  of 
amplitudes  squared. 


2. 5  Types  Of  Noise  and  Their  Distributions 

The  three  main  types  of  noise  distributions  that  will  be  examined  in  this  paper 
are  Poisson,  Gaussian,  and  negative  binomial.  Each  noise  type  is  dependent  upon 
either  the  sensor  conhguration  or  the  source  of  the  target. 

2.5.1  Gaussian.  The  probability  density  function  (pdf)  for  a  continuous 
Gaussian  distributed  random  variable  x  is  [9] 

fx{x)  =  —  00  <  X  <  00,  a  >  0  (6) 

V  271(7 

where  m  is  the  mean  and  a  is  the  variance. 
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Figure  1:  Gaussian  distribution  with  mean  =  0  and  variance  =  1. 


When  the  majority  of  the  noise  is  contributed  by  heat  or  electron  leakage  of  the 
imaging  system,  the  noise  leans  toward  a  Gaussian  distribution.  A  charge-coupled 
device  (GGD)  is  a  device  for  the  movement  of  electrical  charges.  It  is  in  these  GGD 
arrays  of  imaging  systems  where  heat  transfer  and  electron  leakage  is  likely  to  oc¬ 
cur.  In  this  case,  because  the  frequency  spectrum  is  continuous  and  uniform  for  all 
frequency  bands,  the  term  white  noise  is  used.  Additive  white  Gaussian  noise  is  the 
most  common  type  of  noise  model  used  in  electronic  systems  [6]. 

2.5.2  Poisson.  The  probability  mass  function  (pmf)  for  a  discrete  Poisson 
distributed  random  variable  x  is  [9] 


(7) 


where  a  is  the  mean  and  the  variance. 


Poisson  pmf 


Figure  2:  Poisson  distribution  with  mean  and  variance  a  =  10.5. 

The  noise  tends  to  have  a  Poisson  distribution  when  the  dominant  source  of 
noise  comes  from  the  target  when  the  target  is  being  imaged  from  a  passive  system. 
As  a  general  rule,  photons  impinging  on  a  surface  arrive  in  a  random  fashion  after 
being  reflected  from  various  surfaces.  Even  though  the  average  number  of  photons 
arriving  per  second  is  constant,  the  actual  number  for  any  given  second  can  vary 
considerably.  Therefore  this  noise,  often  referred  to  as  shot  noise,  is  modeled  with  a 
Poisson  distribution  [12]. 

2.5.3  Negative  Binomial.  The  pmf  for  a  discrete  negative  binomial  dis¬ 
tributed  random  variable  x  is  [9] 

P[X  =  k]=pk  =  fc  =  r,r  +  l,...  (8) 

where  r  is  a  positive  number  representing  the  number  of  successes  in  the  same  number 
of  independent  Bernoulli  trials  and  p  is  the  probability  of  success.  For  example,  r  =  5 
refers  to  5  successes  in  5  independent  trials. 
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Negative  Binomial  pmf 


Figure  3:  Negative  binomial  distribution  with  probability  p  =  0.4  and  r  represents  10 
independent  Bernoulli  trials. 

Data  collected  from  active  imagers,  such  as  Light  Detection  and  Ranging  (LI- 
DAR)  systems,  appear  to  have  a  negative  binomial  distribution  caused  by  laser  speckle 
resulting  from  the  mutual  interference  of  sets  of  wavefronts  acting  both  constructively 
and  destructively.  The  waves  of  light  from  the  laser  are  coherent  before  they  encounter 
the  target.  Since  most  targets  are  not  perfectly  flat,  the  returning  laser  pulses  are  no 
longer  in  phase.  This  means  that  the  crests  and  troughs  of  the  waves  no  longer  line 
up  with  each  light  pulse.  The  scattered  waves  now  interfere  constructively  to  make 
bright  points  and  destructively  to  make  dark  points.  This  light  and  dark  pattern  is 
referred  to  as  laser  speckle  [2]. 

2.6  Richardson- Lucy  Algorithm 

The  Richardson-Lucy  algorithm  is  an  iterative  deconvolution  method  for  recov¬ 
ering  an  image  that  has  been  blurred  by  a  known  point  spread  function  [7] .  The  image 
estimate  is  updated  upon  each  iteration  of  the  algorithm. 


di  —  ^^PijUj 

j 


(9) 
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The  objective  is  to  determine  the  most  likely  image  estimate,  Uj,  given  the  observed 
image,  di,  and  known  PSF,  pij,  where  u  is  an  estimate  of  u.  This  leads  to  an  equation 
for  Uj  which  can  be  solved  iteratively  according  to 


where 


(10) 


C(  = 

j 

It  has  been  demonstrated  that  if  this  iteration  converges,  it  will  converge  to  the 
maximum  likelihood  solution  [7]. 


2. 7  Related  Work 

While  there  is  much  research  in  the  area  of  blind  deconvolution  of  astronomical 
data,  most  researchers  are  trying  to  solve  the  problem  using  the  short  exposure  atmo¬ 
spheric  model  which  negates  phase  delays  due  to  atmospheric  tilt  [8].  This  research, 
however,  will  examine  the  effects  of  both  the  long  and  short  exposure  optical  trans¬ 
fer  functions  (OTFs).  The  diagram  below  is  a  visual  depiction  of  what  is  going  on 
when  imaging  an  object  through  our  atmosphere.  Equations  in  following  subsections 
explain  the  math  behind  the  diagram. 
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Figure  4:  Astronomical  imaging  diagram 


2.7.1  MAP  Estimation.  The  MAP  estimation  method  in  [4]  uses  an  iter¬ 
ative  optimization  approach  along  with  a  prior  distribution  over  the  quantity  being 
estimated,  to  determine  the  estimate  of  the  scene  that  maximizes  the  likelihood  of  it 
matching  the  true  scene.  The  iterative  approach  of  the  MAP  estimation  algorithm 
researched  for  this  paper  relied  on  the  average  short  exposure  transfer  function,  below, 
to  model  the  atmospheric  effects. 


(12) 


The  average  OTF  of  the  combined  optical  system  and  atmosphere,  Hgys,  used  in  this 
paper,  can  be  modeled  by 


H,y,(x,y)  =  Hyy,(x,y)Hs(x,y) 


(13) 


where  Hopt  is  the  OTF  of  just  the  optical  system  and  Hs  is  the  short  exposure  OTF 
described  by  equation  (12). 
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A  pixel  of  the  noisy  blurred  image  may  be  expressed  as  the  discrete  convolution 
of  Hsys  with  the  remote  scene,  o(a;,  y) 

i{x,y)  =  Hsys{x,y)o{x,y) 

The  probability  of  the  detected  image  pixel  d{x,y)  given  a  particular  observed 
remote  scene  pixel  o{x,y)  can  be  described  by  (14).  The  probabilities  are  determined 
by  the  the  data  distributions  as  seen  in  Section  2.5.  This  particular  probability  has  a 
negative  binomial  distribution. 

P[d{x,y)  =  D{x,y)\o{x,y)  =  0{x,y)] 

{d{x,y)  +  M)\  /  M 

(d(a;,i/)  +  l)!M!V  t{x,y))  V  M  7  ^  ^ 

where  o  is  the  estimate  of  the  target  object  for  a  single  iteration,  M  is  the  coherence 
parameter  of  the  laser  light,  d{x,y)  represents  the  detected  image  and  i{x,y)  is  the 
added  noise  data  [4]. 


MRc) 


(-^0 l^avg) 


(15) 


is  the  assumed  probability  density  function  of  the  a  posteriori  data  based  on  obser¬ 
vation.  is  a  scaling  parameter  and  is  equal  to  the  number  of  pixels  in  the  image. 
Rq  is  the  seeing  parameter  and  Vavg  represents  the  average  value  of  i?o  [4]- 

It  is  from  taking  the  log  of  the  sum  of  equations  (14)  and  (15)  that  the  total 
log  likelihood  function  [4]  is  formed: 


L(o,  i?o) 


f  {d{x,y)  +M)\  \ 
\{d{x,y)  +  l)\M\) 


d{x,y)  X 
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It  can  be  shown  that  the  likelihood  function  [4]  reduces  to 


(16) 


L{o,Ro)  =  d{x,y)ln[i{x,y]  -  [d{x,y)  +  M]  x 


ln[i{x,y)  +  M]j  —N"^— - In 


•  avg 

W 


(17) 


The  maximization  step  [4]  in  which  the  derivative  of  L{o,  Rq)  with  respect  to  o 
is  taken  and  set  to  zero,  yields  the  following 


viv  viv  (d{x,y)\- 


y=\M  +  i{x,y) 


(18) 


An  update  equation  [4]  can  now  be  formed  for  the  seeing  parameter,  similar  to 
the  Richardson-Lucy  approach. 


^^=l^y=l  J  hsysi^  -x,v-y) 

^Li^y=i  (SslJjy)  hsys{x  -^,y-v) 


o^^^{Ro)  =  o°"'(Ro) 


(19) 


where 


i{x,  y)  =  y)  =  -  i,y  -  r/)o°'‘^(C,  y) 

The  MAP  estimation  [4]  now  becomes 


(20) 


Sf=iS^=ihsys(a;  -i,y-y) 


o"""(Ro)  =  o^^\Rq) 


(21) 


2.7.2  APEX.  The  APEX  method  [5]  is  a  general  solution  to  a  specihc 
limited  class  of  blur.  In  order  for  the  technique  to  be  effective,  the  blur  must  be 


14 


symmetric  along  with  certain  other  mathematical  characteristics.  APEX  is  based  on 
a  major  simplifying  assumption  instead  of  an  iterative  approach.  It’s  fast  and  it  does 
not  need  to  know  the  PSF  because  the  algorithm  calculates  it  directly  from  the  image. 
The  equation  for  calculating  the  PSF  directly  from  the  image  is 


h{^,ri)  =  [  (22) 

J 

=  0<a<l,  0</3<l 

where  the  a  and  (5  values  are  determined  by  the  the  specihc  application  of  the  blurring 
function.  For  example,  medical  imagery  and  astronomical  imagery  would  use  different 
values  for  a  and  /3  to  sharpen  the  image.  It  does  require  prior  knowledge  of  these  a  and 
values  used  to  blur  the  image.  For  example,  the  diffraction-limited  optical  transfer 
function  (OTF)  for  a  perfect  lens  can  be  approximated  with  /9  =  |  and  where  a  is 
a  function  of  the  cutoff  frequency.  For  long-exposure  imaging  through  atmospheric 
turbulence,  /^  =  |  and  a  values  are  determined  by  atmospheric  conditions  [8]. 

2.7.3  Results.  The  following  sets  of  images  in  Figures  5  and  6  show  how  the 
MAP  estimation  technique  works  on  a  fully  illuminated  scene  of  a  courtyard.  Fully- 
illuminated  refers  to  the  fact  that  a  majority  of  the  pixels  in  the  image  have  a  value 
greater  than  zero  and  are  discernible  against  a  black  background. 
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(a) 


(b) 


Figure  5:  (a)  is  the  truth  data  and  shows  the  scene  of  the  courtyard  before  any 
distortion,  (b)  shows  the  scene  with  with  blurring  applied.  The  Rq  value  of  the 
simulated  atmosphere  is  10  cm  and  plot  of  the  OTF  can  be  seen  in  the  Figure  6. 


(a) 


(b) 


Figure  6:  (a)  shows  the  side-view  of  the  OTF  that  was  convolved  with  the  courtyard 
image  to  cause  the  blur,  (b)  depicts  the  recovered  scene  from  the  blurred  image  after 
1000  iterations  of  the  MAP  estimator. 


Now  that  this  method  has  been  shown  to  have  satisfactory  results  using  fully  il¬ 
luminated  scenes,  let’s  examine  its  ability  to  accurately  estimate  the  seeing  parameter 
on  astronomical  datasets. 
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(a)  (b) 

Figure  7:  (a)  shows  what  the  binary  star  data  looks  like  without  any  blur.  The  stars 
are  separated  from  each  other  by  1  pixel,  (b)  shows  the  same  binary  stars  blurred 
with  the  same  OTF  in  Figure  6. 
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Figure  8:  This  image  shows  the  difficulty  that  the  MAP  Estimation  algorithm  had 
while  trying  to  recover  the  scene. 

The  difficulty  the  MAP  estimator  has  while  processing  sparse  data  can  be  seen 
in  the  image  above.  It  is  due  to  the  way  the  a  posteriori  information  is  used.  Recall 
the  N^  scaling  factor  from  equation  (15).  It  turns  out  that  this  scaling  factor,  based 
on  the  number  of  pixels  in  the  scene,  is  only  accurate  for  fully  illuminated  scenes. 
In  the  case  of  astronomical  data  sets  where  the  actual  data  is  sparse,  the  scaling 
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factor  drives  the  a  posteriori  data  to  zero,  making  it  useless.  This  can  be  seen  in  the 
following  image. 


Figure  9:  This  image  shows  the  difficulty  that  the  MAP  Estimation  algorithm  had 
while  trying  to  recover  the  scene.  This  plot  shows  that  the  log-likelihood  value  which 
is  way  off  from  the  10  cm  that  it  should  be.  Here  the  weighting  factor,  is  set  to 
the  number  of  pixels  in  the  scene,  16,  384. 

The  next  two  images  in  hgures  10  and  11  show  how  varying  the  value  of  N 
can  result  in  better  estimations.  The  correct  value  for  Rq  is  10  cm  and  is  somewhere 
between  the  estimates  seen  in  these  examples.  The  hrst  guess  of  using  a  scaling  factor 
of  10  does  not  plateau  but  seems  to  increase  asymptotically.  The  next  guess  of  50 
causes  the  estimator  to  converge  prematurely.  One  can  obtain  the  correct  value  of  Rq 
by  messing  with  these  settings  manually  over  a  period  of  trial  and  error.  However, 
these  results  are  image  dependent  and  are  not  a  viable  means  for  automatic  estimation. 


18 


Figure  10:  This  image  shows  result  of  the  MAP  Estimation  algorithm  when  N  was 
set  to  10. 


Figure  11:  This  image  shows  result  of  the  MAP  Estimation  algorithm  when  N  was 
set  to  50. 

The  APEX  method  is  not  without  its  own  difficulties  when  dealing  with  astro¬ 
nomical  data.  Figures  12  and  13  demonstrate  how  varying  the  space  between  the 
stars  has  a  drastic  effect  on  the  ability  of  the  algorithm  to  accurately  estimate  Rq. 
Figures  12  and  13  take  the  row  of  interest  containing  the  binary  stars  and  perform  a 
Fourier  transform  to  examine  the  spatial  frequency  content.  Next  they  compare  the 
amplitude  of  the  spatial  frequencies  to  the  result  of  the  APEX  algorithm. 
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(a)  (b) 


(c) 

Figure  12:  (a)  shows  the  blurry  image  of  a  binary  star  pair  separated  by  one  pixel  (b) 
is  the  result  of  taking  a  slice  of  the  image  through  the  binary  stars  and  plotting  its 
amplitude  (c)  shows  the  result  of  taking  the  Fourier  transform  of  the  data  (subhgure 
b)  and  applying  the  APEX  algorithm.  Both  the  spatial  frequency  of  the  true  data 
and  the  APEX  estimation  are  shown  for  comparison. 
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(a)  (b) 


(c) 

Figure  13:  (a)  shows  the  blurry  image  of  a  binary  star  pair  separated  by  four  pixels 
(b)  is  the  result  of  taking  a  slice  of  the  image  through  the  binary  stars  and  plotting  its 
amplitude  (c)  shows  the  result  of  taking  the  Fourier  transform  of  the  data  (subhgure 
b)  and  applying  the  APEX  algorithm.  Both  the  spatial  frequency  of  the  true  data 
and  the  APEX  estimation  are  shown  for  comparison. 

The  estimation  in  Figure  12  was  accurate  but  it  already  started  to  deviate  from 
the  truth  data  by  the  time  the  stars  were  separated  by  three  pixels,  as  seen  in  Figure 
13.  This  is  due  to  the  fact  that  objects  really  close  to  one  another  appear  as  a  single 
point  source.  As  soon  as  the  separation  becomes  apparent  to  the  APEX  algorithm, 
it  has  difficulty  with  multiple  point  sources. 

MAP  estimation  was  demonstrated  on  fully  illuminated  photographs  and  both 
the  MAP  estimation  and  APEX  blind  deconvolution  methods  were  applied  to  as- 
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tronomical  data  sets.  These  simulations,  both  with  full-scene  and  sparse  datasets, 
looked  at  only  one  particular  optical  system.  The  telescope  receiver  aperture  was 
20  centimeters  and  an  atmospheric  seeing  parameter  of  10  centimeters  was  assumed. 
In  this  particular  case,  the  MAP  estimation  performed  well  on  photographic  images. 
However,  both  methods  were  shown  not  to  be  reliable  when  estimating  astronomical 
data  sets,  each  for  their  own  reasons. 

MAP  estimation  wasn’t  accurate  because  the  scaling  factor  of  the  a  posteriori 
calculation  fails  in  the  absence  of  fully  illuminated  data.  As  soon  as  there  appear  to 
be  more  than  one  point  source,  the  APEX  method  breaks  down  as  well.  It’s  for  these 
reasons  that  the  two  methods  explored  in  this  chapter  are  not  an  effective  means  for 
deblurring  astronomical  data. 

Instead,  Chapter  III  proposes  a  new  technique  in  which  no  information  about 
the  point  spread  function  or  seeing  parameter  is  needed.  This  is  the  usual  case 
when  dealing  with  astronomical  data.  The  proposed  research  effort  would  attempt  to 
develop  a  blind  deconvolution  method  that  produces  useful  results  when  processing 
astronomical  scenes.  The  new  technique  will  compare  the  statistics  of  a  blurry  image 
to  an  estimate  of  the  actual  image,  in  an  iterative  approach,  in  order  to  determine 
the  true  atmospheric  seeing  parameter. 
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III.  Research  Methodology 

This  chapter  deals  with  the  convergence  of  variance  method.  It  briefly  details  the 
approach  below  to  improve  upon  current  work  in  the  area  of  blind  deconvolution. 
Next,  the  expectation  maximization  steps  are  discussed  and  derivations  for  both  the 
Poisson  and  Gaussian  noise  cases  are  provided.  Next  the  coarse  and  hue  search 
approach  to  hnding  the  best  estimated  Rq  is  discussed.  Finally,  the  areas  of  interest 
which  were  covered  under  this  research  for  the  convergence  of  variance  method  are 
discussed.  Example  data  results  for  these  areas  are  provided  in  chapter  4.  These  areas 
include  the  effect  of  SNR,  sparsity  of  the  scene,  varying  seeing  parameter  values,  and 
exposure  type. 

3.1  Convergence  of  Variance  Method 

The  convergence  of  variance  (COV)  method  that  was  developed  for  this  research 
has  two  main  parts.  The  hrst  part  of  this  method  involves  deriving  the  expectation 
maximization  equations  needed  for  a  particular  noise  distribution  of  the  data.  This 
will  allow  for  the  correct  estimations  of  the  images  to  be  made.  The  next  part  of  this 
method  involves  using  the  variance  of  the  observed  image  as  the  stopping  criteria  for 
the  expectation  maximization  algorithm. 

The  hrst  step  is  to  acquire  an  astronomical  image  and  compute  its  variance. 
Next,  in  an  iterative  Richardson-Lucy  approach,  the  variance  of  the  observed  image 
and  the  variance  of  the  latest  estimate  of  the  image  are  compared.  However,  unlike 
in  the  Richardson-Lucy  approach,  here  the  PSF  does  not  need  to  be  known.  For 
the  hrst  estimate,  the  seeing  parameter  estimate  is  set  equal  to  1.  This  process  will 
continue  until  the  length  of  the  iterations,  needed  to  converge,  reaches  a  predetermined 
threshold  set  by  the  image  analyst.  Once,  this  threshold  is  met,  the  value  of  the  seeing 
parameter  estimate  is  incremented  and  the  process  repeats.  Once  the  variance  of  the 
image  estimate  converges  with  the  variance  of  the  acquired  astronomical  image,  the 
best  estimate  of  the  seeing  parameter  is  known  and  the  image  can  be  seen  as  it  would 
appear  without  the  blurring  ehects  of  the  atmosphere. 
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3.2  Expectation  Maximization  (EM)  Algorithm 

An  EM  algorithm  is  widely  used  in  the  iterative  calculations  of  maximum  likeli¬ 
hood  estimates.  The  EM  algorithm  is  not  a  single  algorithm  but  rather  a  framework 
or  approach  to  solving  a  problem  [3].  During  each  iteration  or  estimate  rehnement, 
the  EM  algorithm  goes  through  both  an  expectation  and  a  maximization  step.  The 
initial  values  or  parameters  are  estimated  iteratively  until  convergence.  Other  good 
sources  for  papers  which  use  an  EM  approach  to  blind  deconvolution  can  be  found 
in  [1,4,10,11].  The  EM  derivations  for  images  containing  Gaussian  and  Poisson 
random  variables  are  included  below.  An  example  of  EM,  using  negative-binomial 
random  variables,  is  found  in  Chapter  2  Section  7.1,  MAP  estimation. 


3.2.1  EM  Steps.  Here  are  the  main  steps  in  the  EM  process.  These  steps 
are  outlined  below  in  each  of  the  derivations. 


1.  Obtain  statistical  model  for  the  measured  (incomplete)  data,  d 

2.  Obtain  statistical  model  for  the  complete  data,  d 

3.  Formulate  the  complete  data  log-likelihood,  L 

4.  Compute  expectation  step 

5.  Compute  maximization  step 

6.  Compute  update  equation 


3.2.2  Poisson  EM  Derivation. 

Step  1:  Obtain  statistical  model  for  the  measured  data: 


E  [d{x,y)] 
P[d{x,y)] 


N 

i{x,y)  =  o{z,  w)h{x  —  z,w  —  y) 

Z,W  =  1 

i{x, 

d(x,y)\ 


(23) 


The  measured  data  is  an  image  in  which  each  pixel  is  a  realization  of  a  Poisson  ran¬ 
dom  variable  whose  mean  is  the  result  of  a  convolution,  h  is  a  linear  operator  that 
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transforms  the  true  image,  o,  into  the  mean  of  the  noisy  blurred  image,  i. 

Step  2;  The  complete  data:  Invent  a  set  of  mythical  (complete)  data,  d  and  a  re¬ 
lationship  between  this  set  and  the  measured  (incomplete)  data,  d.  We  choose  the 
complete  data  to  be  related  to  the  incomplete  data  through  equation  (24). 

d{x,y)  =  ^d{x,y,z,w)  (24) 

Z,W 

Statistical  model  for  the  complete  data:  Select  a  statistical  model  for  the  complete 
data  such  that  it  produces  the  statistical  model  for  the  incomplete  data  through  their 
dehned  relationship.  We  chose  the  complete  data  to  be  Poisson  as  well  since  the  sum 
of  Poisson  random  variables  is  also  Poisson. 

E  ^d{x,y,  z,w)  =  o{z,w)h{x  —  z,y  —  w) 

E  ^d{x,y,z,w) 

_  Z,W 

^o{z,w)h{x  -  z,y  -  w)  =  E[d{x,y)]  (25) 

Z,W 
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Step  3:  Formulate  the  complete  data  log-likelihood: 


P  d{x,y,  z,w) 

d 


L(o,  h) 


o{z,  w)h{x  -  z,y- 

d{x,y,z,w)\ 

d{x,  y,  z,  w)\f{x,  y,  z,w)  e  I (1,  N) 
o{z,  w)h{x  -z,y- 


x,y,z,w=l 


d{x,y,z,w)\ 


In  (^P  ^d{x,y,z,w)\/{x,y,z,w)  e  I{l,N)  ) 

N 

d{x,y,  z,w)\n{o{z,w)h{x  —  z,y  —  w)) 


x,y,z,w=l 


o{z,  w)h{x  —  z,y  —  w)  —  \n{d{x,  y,  z,  w)\) 


Step  4:  Find  the  expected  value  of  the  complete  data  log-likelihood. 

E  [L{o,h)\d{x,y),o°^‘^,h]  =  Q{o,h) 

=  ^  E  d{x,y,z,w)\d{x,y),o°^‘^,h 

x,y^z,w 

*  [ln(o(;2,  w)h{x  —  z,y  —  tc))]  —  o{z,  w)h{x  —  z,y  —  w) 

—  E  \n{d{x,y,  z,w)\\d{x,y),o°’‘‘^  (27) 

Dehne  two  statistically  independent  Poisson  variables  that  adhere  to  the  relationship 
dehned  in  (24). 

di  =  d{x,y,Zo,Wo)d2 

=  '^d{x,y,z,w)  -  d{x,y,zo,wo)  (28) 

Z,W 

Dehne  the  expected  value  and  probability  of  Poisson  randomly  distributed  datasets. 
Here  is  the  simple  case  of  two  arbitrary  Poisson  variables.  It’s  easier  to  determine 
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the  probability  when  summing  two  Poisson  variables,  than  it  would  be  for  a  Poisson 
variable  for  each  pixel  in  the  scene.  Since  the  sum  of  Poissons  is  also  Poisson,  this 
case  can  be  extented  to  n  Poisson  variables. 


E[d^] 

E[d2] 

d{x,y) 

P[di,d‘^ 

(^2 

P[(i,  di] 


mi 


m2 


d  —  di  (^2 


di\ 
d  —  di 


do! 


di,  d2  =  0, 1,  2, 


m^  e 


di!  (d-di)! 


(29) 


Next,  hnd  the  probability  of  the  complete  data  given  the  measured  dataset.  This 
is  the  case  in  which  each  pixel  is  a  Poisson  random  variable  and  d  is  the  sum  of  all 
variables. 


F[di|d]  =  P[d,di]/P[d] 


=  m 


P[di|d] 

P[di|d] 


^{d  di)^- 


m2 


P[di|d]  = 


di! 

(d!)mf^ 


(d-di)! 


(mi  +  m2Ydi\  (d  —  di 


(d!)m^ 


di 


m 


{d-d-i) 


(mi  +  m2)‘!'^i+'^“‘^ddi!  (d  —  di)! 

(d!)(mi/[mi  +  m2])'^^  (m2/[mi  +  m2])‘''^“'^^^ 


di! 

d(mi/[mi  +  m2]) 


(d-di)! 


(30) 


E 


d{x,y,z,w)\d{x,y),o°  ,h 


d{x,  y)o°^'^{z,  w)h{x  —  z,y  —  w) 


4 old  I 


x,y) 


(31) 
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Now  substitute  (31)  into  Q(o,h)  from  (27). 


Q{o,h)  =  ^ 


d{x^  w)h{x  —  z^y  —  w) 


x,y,z,w 


iold(^X^  y) 

*  [ln(o(2:,  w)h{x  —  z,y  —  w))]  —  o{z,  w)h{x  —  z,y  —  w) 

-  E[\n{d{x,y,z,w)\)\d{x,y,oy’^'^,h]  1  (32) 


Step  5:  Maximize  the  result  by  first  taking  the  derivative  and  then  setting  it  equal  to 
zero.  The  variables  Zq  and  Wq  are  specihc  locations  of  the  generic  variables,  and  w. 


dQ{o,  h) 
do{zo,wo) 


E 


x,y,z,w 


d 


do{zo,wo) 


d{x,y)o°^^{z,w)h{x-z,y-w) 

*  ^ ^ - m[o{z,w)h[x  —  z,y  —  w)) 


d 


do{zo,wo) 

d 

do{zo,WQ) 


iold(^X^  y) 
o{z,  w)h{x  —  z,y  —  w) 

E[\n{d{x,  y,  z,  w)\)\d{x,  y)o°^‘^] 


(33) 


Now  that  the  derivative  has  been  taken,  most  terms  go  away  and  the  remaining  term 
is  set  equal  to  zero 


0  = 


E 


x,y,z,w 


d 


do{zo,wo) 


E[ln{d{x,y,z,w)\)\d{x,y)o' 


oldi 


d 


do{zo,wo) 


o{z,w)h{x  —  z,y  —  w)  =  S{z  —  zo,w  —  wo)h{x  —  z,y  —  w) 


d 

do{zo,wo) 


* 

* 


d{x,  w)h{x  —  z,y  —  w) 

y) 

[111(0(2;,  w)h{x  —  z,y  —  w))] 
d{x,  y)o°’''^{z,  w)h{x  —  z,y  —  w) 
iold(^X^  y) 

d{x,  y)o°^‘^{z,  w)h{x  —  z,y  —  w) 

iold(^X^  y) 

d{x,  y)o°^'^{z,  w)h{x  —  z,y  —  w) 

iold(^X^  y) 


1  do{z,w) 
do(z,  w)  do{zo,  Wo) 

- r(5(2:  -  2:0,  W  -  Wo) 

ao(2:,  w) 


(34) 


Maximize  Q  with  respect  to  h{zo,Wo). 


dQ{o,h) 

dh{zo,Wo) 


d{x^  w)h{x  —  z^y  —  w) 


rioldl 


X,  y) 

x,y,z^w  '■ 

5{x-z-Zo,y-w-  Wq) 
h{x  —  z,y  —  w) 


S{x  —  z  —  zo,y  —  w  —  wo)o{z,  w) 

(35) 


Now,  apply  the  sifting  property  of  the  Dirac  function: 


dQ{o,h) 

dh{zo,Wo) 

dQ{o,h) 

dh{zo,Wo) 


0 


E 


d{x,y)o°’‘^{z,w)h{x  —  z,y  —  w)  1 

i^^<^{x,y)  h{zo,Wo) 


o{x  -  zo,y-wo) 


0 


E 


d{x,y)o°’‘^{z,w)h{x  —  z,y  —  w)  1 

i^^<^{x,y)  h{zo,Wo) 


o{x  -  Zo,y-Wo) 


(36) 
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Find  the  new  estimate  for  the  object 


^o{x  -  Zo,y-  Wo) 
x,y 


h{zo,Wo) 


o^^^{zo,wo) 


E 


d{x,  y)o°^‘^{zo,  wo)h{x  —  z,y  —  w) 


jold^X^  y) 


h{zo,wo) 


h(-y  in  \  d{x,y)o°^‘‘{z,w)h{x-zo,y-wo) 

IT'KZO:  Wo)  y  iold(^^^y) 

Y.x,yO{.x-zo,y-wo) 

^oldf  ^  ^  d{x,y)h{x-zo,y-wo) 

O  [Zo,  Wo)  l^^  y - iSny^) - 

E.,j  l‘{x-  zo,y  -  Wo) 


(37) 


Step  6:  Form  the  update  equation  for  estimated  object 

We  chose  h  summing  to  1  to  be  a  constraint  since  this  will  ensure  the  conservation  of 
energy  in  the  system.  That  is,  both  the  true  image  and  the  estimated  image  would 
both  contain  the  same  number  of  photons. 


Choose  h{x  —  zo,y  —  wo)  =  1 


x,y 


Show  that  the  OTF  update  sums  to  1 


(38) 


h(zo,wo)  = 


J^h(zo, 


Wo  = 


x,y 


h(z  in  d(x,y)o°‘'^(x-zo,y-WQ) 

'^[Xo,Wo)  2_^^  y  iold(^x,y) 

T.x,yOi^  -  -  Wo) 

h(z  in  \  d(x,y)o°‘‘^{x-zo,y-wo) 

'^[Xo,Wo)  l^^  y  iold(x^y) 

Y.x,yd{x,y) 

ST  hf-y  in  d{x,y)o°^‘^{x-zo,y-wo) 

l^x,y  WO)  Z^x,y  i°‘d)x,y) 

J2x,yd{x,y) 

T.x,y  T.x,y  -Z0,y-  Wo)h{zo,  Wo) 


J2x,yd{x,y) 


(39) 
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(40) 


J2h^^^{zo,Wo) 

x,y 


^  -  Zo,y-  Wo)h{zo,  Wq) 

x,y 

T,x,yd{x,y)  ^  ^ 

J2x,yd{x,y) 


3.2.3  Gaussian  EM  Derivation. 

Step  1:  Obtain  statistical  model  for  the  measured  data: 


E  [d{x,y)] 
P[d{x,y)\ 


N 

i{x,y)=  ^  o{z,w)h{x  -  z,w  -  y) 

Z,W=1 

1  (d(x,y)-i(x,y))^ 

, _  e  ^ 

y/^a 


(41) 


The  measured  data  is  an  image  in  which  each  pixel  is  a  realization  of  a  Gaussian 
random  variable  whose  mean  is  the  result  of  a  convolution,  h  is  a  linear  operator  that 
transforms  the  true  image,  o,  into  the  mean  of  the  noisy  blurred  image,  i. 


Step  2:  The  complete  data:  Invent  a  set  of  mythical  (complete)  data,  d  and  a  re¬ 
lationship  between  this  set  and  the  measured  (incomplete)  data,  d.  We  choose  the 
complete  data  to  be  related  to  the  incomplete  data  through  equation  (41). 


N 

d{x,y)  =  dix,y,z,w)  (42) 

z,w=l 

Statistical  model  for  the  complete  data:  Select  a  statistical  model  for  the  complete 
data  such  that  it  produces  the  statistical  model  for  the  incomplete  data  through  their 
dehned  relationship.  We  chose  the  complete  data  to  be  Gaussian  as  well  since  the 
sum  of  Gaussian  random  variables  is  also  Gaussian. 
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d{x,y,z,w) 
P  d{x,y) 


=  o{z,w)h{x  —  z,w  —  y) 

X  {d{x,y,z,w)  —  o{z,w)h{x  —  z,w  —  y))^ 


Step  3:  Formulate  the  complete  data  log-likelihood: 


d  =  d{x,y,  z,w)\/{x,y,  z,w)  E  N) 

N  . 

~  "1  -m — w-  \_  {d{x,y,z,w)  —  o{z,w)h{x  —  z,w  —  y)) 

P  d  =  n 

x^y,z,w=l  ^ 

L{o,h)  =  In  d{x,y,  z,w)W{x,y,  z,w)  E  I{1,  N) 


{d{x,  y,  z,  w)  —  o{z,  w)h{x  —  z,w  —  y)Y 


x^y,z^w=l 


Step  4:  Find  the  expected  value  of  the  log-likelihood. 


Q{o)  =  E[L{o)\d{x,y),Ooid{z,w)]  =  ^  y- E[{d{x,y,  z,w)\ooid{z,w),d{x,y)f] 

x,y,z,w 

+  2E[{d{x,  y,  z,  w) \ooid{z,  w),  d{x,  y))] 

+  o{z,w)'^h{x  —  z,y  —  (45) 

Now  hnd  this  piece  of  (46)  by  dehning  two  sets  of  Gaussian  random  variables  in  (48). 


d{x,y,z,w)\d{x,y),o{z,w) 


First,  the  two  sets  of  Gaussian  randomly  distributed  data  are  dehned  below. 


—  di  d2 
=  d{x,y,zo,wo) 

N 

=  ^  d{x,  y,  z,  w)  -  d{x,  y,  zq,  wq) 
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The  mean  and  variances  for  the  complete  data  are  provided  below. 


(48) 


^[di] 

=  mi,  V ar[di]  = 

P[d2] 

=  m2,  V^ar[d2]  =  a‘^{rd  —  1) 

d{x,y) 

=  d  =  di  “h  g?2 

(49) 

Now  the  probability  of  complete  data  is  calculated 


P[di,  d2 


P[d^  di 


1  (  1  (  {d2-m2f\ 

spHo  V  2^2  J  a^27r(n^  -  I)""'*’  V  2<T^(n^  -  1)  j 

1  /  (di-mi)^\  (  {d-di-m2y 

I - ^  j  “P  I" 


(50) 


The  probability  of  complete  data  can  be  determined,  given  the  set  of  measured  data 


P[di|d] 

P[di\d] 

P[di\d] 

P[di\d] 

E[di\d] 


27rfT^\/n^  — 1 


exp 


2^2 


exp 


(d— di— m2)2 
2(72(n2  — 1) 


\/2'ko 


72^ 


n 


-exp 


exp 


*  exp 


2a\/n{n?  —  1) 

(di  -  mif 


2u2 


n 


2(T\J'x{'n?  —  1) 

y/2n 


exp 


{d—di-m2y 

2<j‘^'nP‘ 


{d  —  di  —  m2Y 


exp 


2a^n‘^ 

{d  —  di  —  7712)^ 


-.exp 


2(TA/vr(7p^^T) 

{d  —  mi  —  m2) 


2a‘^{rP  —  1) 

{in?di  —  n^mi  +  m2  +  rui  — 

2'n?a‘^['n?  —  1) 

^^4 


=  mi  + 


m 


(51) 
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Substituting  these  results  back  into  the  previous  form  of  the  equation  yields  the 
following  result. 


E[d{x,  y,  z,  w)\ooid{z,  w),d{x,  y)] 


Ooid{z,w),  h{x  -  z,y  -w) 

-{ioid{x,  y)  -  Ooidiz,  w)h{x  -  z,y  -  w)) 


Ooid{z,  w)h{x  -  z,y-w)  +  d{x,  y) 


Ooid{z,  w)h{x  -  z,y  -w)  + 


{d{x,y)ioid{x,y)) 


(52) 


Step  5:  Maximize  the  result  by  taking  the  derivative  and  setting  it  equal  to  zero. 


dQ{o) 


E 


dQ{o) 


do{zo,wo) 


-  E[{d{x,  y,  z,  w)\ooid{z,  w),d{x,  y)f 


+  2E[{d{x,  y,  z,  w)\ooid{z,  w),  d{x,  y))] 
+  o{z,wYh{x  —  z,y  —  wY 


o{z,  w)h{x  —  z,y  —  w) 


(53) 


Equation  (55)  is  the  first  piece  of  (54)  and  doesn’t  depend  on  o{zo,Wo),  so  taking 
and  setting  it  equal  to  0  yields: 

ao{ZQ,WQ)  on  J 


E 


x,y^z^w 


d 


do{zo,wo) 


E[{d{x,y,z,w)\ooid{z,w),d{x,y)f]  =  0 


(54) 


Equation  (56)  is  the  second  piece  of  (54). 


E 


x^y,z,w 


d 


do{zo,  Wo) 


2o{z,  w)h{x  —  z,y  —  w) 


Ooid{z,  w)h{x  -  z,y  -w)  + 


{d{x,y)  -i{x,y)) 
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E 

x,y,z,w 


2o{z 


zq,  w  —  wo)h{x  —  z,y  —  w)* 


(^Ooid{z,  w)h{x  -  z,y  -w)  + 


{d{x,y)  -i{x,y))\ 
ri^  ) 


E 


x,y  L 


2ooid{z,  w)h{x  -  zo,y  -  wof  + 


h{x  -  Zo,y-  Wo)  * 


{d{x,y)  -i{x,y)) 


(55) 


Equation  (57)  is  the  third  and  hnal  piece  of  (54). 


E 


x,y^z,w 


do{z,  w)‘^h{x  —  z,y  —  w)" 
do{zo,wo) 


\^d{z  —  zo,w  —  wo)o{z,  w)h{x  —  z,y  —  wY 
^  [20(2:0,  Wo)h{x  -  Zo,y  -  WoY 


x,y,z,w 


Substituting  the  results  of  (55), (56),  and  (57)back  into  (54)  yields: 


(56) 


0  = 


dQ{o) 

do{zo,wo) 


E 


2ooid{zo,  wo)h{x  -  zo,y  -  wq)"^ 


2n2 


h{x  -  zo,y-  Wo) 


(d(x,y)-i(x,y)) 


2a^ 

/ 20(2:0,  wo)h{x  -  zo,y  -  wo)^ 

V  2a2 


E 


20(2:0,  wo)h{x  -  zo,y-  WoY 


2^2 

^,y 

h{x  -  zo,y  -  Wo) 


E 


x,y 


2ooid(zo,  wo)h(x  -  zo,y-  wo)" 
2n2 


(d(x,y)-i(x,y)) 

r)  2 


2a2 


+ 


(57) 
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Step  6:  This  results  in  the  following  update  equation: 


Onew{z,w)  =  Oold{z,w)  + 


E 


N 

X  =  1 


T.y=iHx  -  Z,y  -  w){d{x,y)  -ioid{x,y)) 

Ef=i  EjLi  hix  -z,y-  w)^ 


(58) 


3.3  Using  Coarse  and  Fine  Convergence 

Using  coarse  convergence  followed  by  hne  convergence  is  a  method  for  quickly 
and  accurately  identifying  the  seeing  parameter  and  estimate  of  the  actual  image. 
During  the  coarse  convergence  step,  the  selected  EM  algorithm  performs  iterative 
image  estimations  while  searching  through  integer  values  or  multiples  of  i?o,  if  the 
user  chooses,  until  the  algorithm  converges.  For  example,  one  could  do  a  coarse  search 
for  Rq  by  multiples  of  two  to  hnd  an  estimate  quicker.  Next,  the  hne  convergence 
begins  by  starting  its  search  just  before  the  coarse  estimation  value  of  Rq.  This  time, 
however,  the  algorithm  performs  iterations  using  decimal  values  of  Rq  to  narrow  in 
on  the  best  estimate. 

This  two  step  approach  serves  two  distinct  purposes.  The  hrst  is  to  quickly 
narrow  down  the  search  area  for  Rq.  The  second  purpose  is  to  rehne  the  search  for 
a  more  accurate  result  for  an  estimated  value  of  Rq,  since  not  all  Rq  values  will  be 
integers.  Figure  14  demonstrates  the  process  of  hrst  searching  through  integer  values 
of  Rq  before  rehning  the  search  area. 


Coarse  Search 


Figure  14:  Diagram  of  the  coarse  and  hne  search  method 
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Figures  15  through  17  demonstrate  this  process  using  simulated  starfield  data. 


True  Object  Blurred/Noisy  Image 


Figure  15:  Short  exposure  blur  caused  by  atmospheric  distortion  with  seeing  param¬ 
eter  of  4.3  cm  and  Poisson  noise  distribution. 

Figures  16  and  17  each  contain  two  images.  The  plot  on  the  left  is  a  snapshot 
of  the  convergence  process  and  the  image  on  the  right  reflects  the  updated  estimate 
of  the  image  based  on  the  level  of  convergence.  The  variance  of  the  estimated  image 
converges  with  the  variance  of  the  measured  image  with  each  iteration  and  the  image 
on  the  right  is  updated.  There  is  also  some  key  data  provided  with  the  convergence 
plots.  The  current  Rq  search  value,  displayed  above  the  plot,  the  mean  squared  error 
and  the  estimated  number  of  iterations  needed  to  converge  are  also  updated  upon 
each  iteration.  If  the  iterations  needed  exceeds  the  user  dehned  level,  the  Rq  search 
value  is  increased  and  the  process  begins  again.  Once  the  iterations  needed  reaches 
less  than  1  the  process  is  complete  and  the  COV  method  displays  both  the  coarse  and 
fine  image  results. 
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Number  Of  Iterations  ( estimation  =  0.17324  ) 


Object  Estimate  Using  Blind  Deconvolution 


X  Coordinate 


Figure  16:  Coarse  search  by  Ro  increments  of  1  passed  up  the  correct  value  of  Ro  = 
4.3  cm  and  instead  stops  at  the  best  coarse  value  of  5  cm. 


Object  Estimate  Using  Blind  Deconvolution 


X  Coordinate 


Figure  17:  Fine  search  by  Rq  increments  of  .1  allows  for  the  correct  value  of  Rq  =  4.3 
cm  to  be  found. 

Here  are  some  more  results  of  the  coarse  and  hne  search  method.  The  image 
below  in  Figure  18  depicts  the  true  scene  along  with  the  image  containing  atmospheric 
distortion,  as  viewed  by  an  optical  system.  Note  that  the  two  stars  in  the  middle 
appear  as  one  in  the  blurred  image. 
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True  Object  Distorted  Image 


Figure  18:  Short  exposure  blur  caused  by  atmospheric  distortion  with  seeing  param¬ 
eter  of  5  cm. 

Figures  19,  20,  and  21  below  show  the  results  of  the  coarse  and  hue  conver¬ 
gence  estimates  on  three  sets  of  Gaussian  distributed  images  with  differing  seeing 
parameters. 
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Figure  19:  Side-by-side  comparison  using  simulated  star  held  data.  Actual  Rq  value 
is  5  cm. 
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Coarse  Search  Fine  Search 


Figure  20:  Side-by-side  comparison  using  simulated  star  field  data.  Actual  Ro  value 
is  10  cm. 


Coarse  Search  Fine  Search 


Figure  21:  Side-by-side  comparison  using  simulated  star  field  data.  Actual  Rq  value 
is  15  cm. 

In  the  first  two  figures,  the  coarse  search  was  right  on  but  in  the  third  it  overshot 
the  correct  value  of  Rq.  The  fine  search  was  able  to  estimate  the  correct  value  to  within 
0.1  cm.  These  variations  are  all  due  to  the  the  SNR  ratios  of  the  images.  Both  the 
coarse  and  fine  search  estimates  are  more  accurate  as  the  SNR  improves. 

3.4  Image  Parameters  Of  Interest 

Some  key  parameters  which  may  effect  the  results  of  the  COV  method  include 
the  SNR  and  sparsity  of  the  data,  the  atmospheric  seeing  parameter  value  at  the  time 
the  image  was  acquired  and  the  exposure  type.  These  parameters  will  be  explored 
more  deeply  along  with  sample  data  to  support  the  ideas  in  Chapter  4.  The  SNR 
of  the  data  may  affect  the  ability  of  the  COV  method  to  accurately  estimate  Rq. 
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Different  values  of  Rq  may  also  bias  the  results.  Sparse  and  full  scene  data  refer  to 
the  amount  of  illuminated  pixels  in  the  scene.  If  the  COV  method  can  produce  valid 
results  for  sparse  data  as  well  as  full  scene  data,  then  it  will  have  shown  improvement 
over  the  MAP  estimation  and  APEX  techniques  explored  in  Chapter  2.  Finally  the 
exposure  time  of  the  images  will  be  examined  to  see  if  it  also  has  any  bearing  on  the 
ability  of  the  COV  method  to  converge  and  determine  the  best  estimate  of  Rq. 
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IV.  Results  and  Analysis 

Chapter  4  examines  the  effects  of  SNR,  sparsity  of  the  scenes,  varying  seeing  parameter 
values,  and  exposure  type  and  provides  example  simulated  data.  It  also  summarizes 
the  results  of  the  COV  method  as  compared  to  the  previous  work  researched  in  the 
MAP  estimator  and  APEX  algorithm.  Next  the  results  are  provided  for  the  COV 
method  used  on  laboratory  measured  data  and  the  idea  of  using  the  COV  method 
on  other  data  types  is  advocated.  Finally,  limiting  factors  in  the  effectiveness  of  the 
COV  method  are  explored. 


4-1  SNR  Effects 

SNR  plays  a  key  role  on  the  convergence  of  the  COV  method.  If  the  SNR  is 
too  low,  the  COV  method  tends  to  converge  relatively  quickly  at  an  Rq  value  which 
is  less  than  the  known  seeing  parameter  of  the  simulated  atmosphere.  This  quick 
convergence  with  few  iterations  at  that  Rq  value  leads  to  heavily  blurred  results. 
As  long  as  the  SNR  is  high  enough,  the  COV  method  does  not  have  an  issue  with 
convergence.  Figure  22  contains  an  image  of  a  courtyard  along  the  same  image  as  it 
would  appear  if  blurred  by  the  short  exposure  OTF  and  with  added  Poisson  noise. 
The  same  effect  can  be  demonstrated  with  the  starfield  datasets,  but  the  lack  of 
illuminated  pixels  in  the  scene  make  it  more  difficult  to  detect  the  differences  by  eye. 


True  Image  Blurred  Image 


Noisy  Image 


(a)  (b)  (c) 

Figure  22:  (a)  shows  the  image  of  the  true  scene,  (b)  contains  the  same  scene  blurred 
with  the  short  exposure  OTF,  and  (c)  shows  the  blurred  image  with  Poisson  noise. 
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Figure  23  shows  the  difference  in  performance  of  the  COV  method  when  the 
SNR  of  the  observed  image  is  varied.  The  SNR  is  scene  dependent  and  the  minimum 
SNR  value  needed  to  converge  properly  was  not  explored  for  the  purposes  of  this 
research.  The  Rq  values  in  Figure  23  are  the  estimated  values.  The  true  value  for  Rq 
is  5  cm. 


Coarse  Search  Coarse  Search 


(a)  SNR  =  5.7  dB  (b)  SNR  =  11.2  dB 

Coarse  Search 


(c)  SNR  =  19.4  dB 

Figure  23:  (a)  The  COV  method  estimated  a  lower  Rq  value  and  converged  early 
without  enough  iterations  to  deblur  the  image  because  the  SNR  was  only  5.7  dB,  (b) 
the  COV  method  hnds  correct  Rq  value  but  converges  after  few  iterations  because  of 
low  SNR  of  11.2  dB,  (c)  SNR  was  increased  to  19.4  dB  which  resulted  in  an  image 
estimate  that  closely  resembles  the  true  image 


4-2  Effect  Of  Different  Seeing  Parameters 

Varying  the  seeing  parameter  values  of  the  simulated  atmosphere  seemed  to 
have  no  real  effect  on  the  results  of  the  image  reconstruction.  The  only  thing  that 
changed  was  the  COV  method  took  a  longer  amount  of  time  to  converge.  The  images 
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below  in  Figure  24  show  the  true  scene  and  the  effects  of  atmospheric  turbulence. 
The  next  set  of  images  in  Figure  25  show  the  results  of  the  same  simulated  star  helds 
when  using  three  different  Rq  values. 


True  Object  Distorted  Image 


Figure  24:  Here  is  the  image  of  the  true  scene  with  no  blurring  and  the  same  scene 
after  being  convolved  with  the  short  exposure  OTF  with  Rq  =  10  cm. 
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(c)  Actual  Rq  =  15cm 

Figure  25:  These  three  images  have  a  Gaussian  distribution  and  were  blurred  with 
a  short  exposure  OTF.  The  SNR  was  held  constant  at  18.6  dB  for  all  three  images. 
The  images  show  the  estimated  Rq  values  after  being  deblurred,  while  the  labels  show 
the  actual  Rq  values. 

4.3  Short  y^s  Long  Exposure 

Exposure  refers  to  the  amount  of  light  gathered  or  the  integration  period  of  the 
image.  Short  exposure  images  are  typically  formed  by  taking  snapshots  of  a  target 
with  less  than  a  second  of  integration  time  for  the  light  gathering  of  the  optical  system. 
Long  exposure  images,  on  the  other  hand,  are  usually  formed  by  collecting  light  for 
many  seconds. 

Short  exposnre  images  of  astronomical  scenes  are  far  clearer  as  the  change  in 
the  seeing  parameter  is  negligible  in  such  a  short  time.  They  also  do  not  contain  any 
blurring  effects  caused  by  the  motion  of  the  earth  and  the  target  in  qnestion. 
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The  image  below  in  Figure  26  depicts  the  true  scene  along  with  the  image, 
containing  atmospheric  distortion,  as  viewed  by  an  optical  system.  Note  that  the  two 
stars  in  the  middle  appear  as  one  in  the  blurred  image. 


True  Object  Distorted  Image 


Figure  26:  Long  exposure  blur  caused  by  atmospheric  distortion  with  seeing  parameter 
of  i?o  =  5  cm. 

The  simulated  data  below  in  Figures  27  and  28  contains  both  short  and  long 
exposure  image  estimates  of  the  same  star  held  data  set.  The  COV  method  attempts 
to  deconvolve  both  sets  of  images  and  estimate  the  true  scene.  The  image  has  a 
Gaussian  distribution  and  the  actual  i?o  value  is  5  cm. 
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Figure  27:  Short  exposure(L),  long  exposure(R).  Side-by-side  comparison  of  the  de- 
blurred  results  using  simulated  star  held  data. 
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Fine  Search 
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Figure  28:  Short  exposure(L),  long  exposure(R).  Side-by-side  comparison  of  the  de¬ 
barred  results  using  simulated  star  held  data. 

The  COV  method  was  capable  to  accurately  estimate  the  short  exposure  image 
with  the  binary  pair  clearly  visible.  The  long  exposure  image,  containing  motion  had 
too  much  blur  and  the  COV  method  was  unable  to  distinguish  between  the  two  stars 
in  this  case.  However,  it  is  brighter  than  the  surrounding  stars  of  the  same  intensity, 
indicative  of  its  larger  size. 

4-4  Sparse  Vs  Full  Scene  Data 

Sparse  and  full  scene  data  refer  to  the  amount  of  illuminated  pixels  in  the  scene. 
For  the  purpose  of  this  research,  an  image  with  over  half  of  its  pixels  illuminated  is 
referred  to  as  full  scene  and  any  image  containing  less  is  referred  to  as  sparse  data. 
Though  there  are  many  images  that  range  in  sparsity,  this  research  focused  on  the 
extremes.  The  data  includes  binary  and  small  star  clusters  as  well  as  photographic 
images  of  a  courtyard. 

Chapter  2  showed  how  the  MAP  estimator  was  effective  on  full  scene  data 
but  and  both  the  MAP  estimator  and  APEX  algorithms  were  not  accurate  when 
processing  sparse  datasets.  On  the  other  hand,  the  COV  method  in  this  research  has 
demonstrated  effectiveness  with  both  full  scene  data  as  seen  in  section  4.1  and  sparse 
data  in  sections  4.2  through  4.3. 
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4-5  Results  of  Laboratory  Data 

Sections  1  through  3  conhrmed  the  ability  of  the  COV  algorithm  where  Chapter 
2  showed  how  the  MAP  estimator  and  APEX  algorithm  failed.  It  featured  simulated 
data  to  produce  results  which  demonstrated  the  desired  effects  for  each  section.  Now  a 
set  of  measured  data  obtained  under  laboratory  conditions  is  provided  to  give  support 
that  the  COV  method  will  work  on  real  measured  data. 

The  data  in  Figure  29  is  intended  to  represent  a  large  celestial  object  such  as 
the  moon.  It  was  captured  by  mounting  an  LED  behind  a  bar  chart  in  a  dark  room. 
The  focusing  lens  was  adjusted  to  mimic  the  effects  of  the  atmosphere  of  the  local 
area  in  Dayton  Ohio,  which  is  typically  between  4  and  6  cm  on  a  clear  night. 


Figure  29:  Target  object  as  observed  by  optical  system 

Since  this  data  is  measured  there  is  no  way  to  adjust  the  SNR,  sparsity  of  the 
scene  or  the  OTF.  The  only  parameter  that  can  be  manipulated  is  the  amount  of 
iterations  that  the  image  analyst  will  permit  the  algorithm  to  perform.  Figures  30 
through  32  show  the  effect  that  the  number  of  iterations  has  on  the  ability  of  the 
algorithm  to  hud  the  best  estimate  of  Rq  and  to  deblur  the  image. 
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Coarse  Search  Fine  Search 


Figure  30:  The  algorithm  converged  late  because  the  maximum  iteration  allowed  was 
too  low.  This  image  pair  demonstrates  the  result  if  the  user  sets  the  max  iteration 
level  to  200. 
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Figure  31:  The  algorithm  still  converged  late  because  the  maximum  iteration  allowed 
was  too  low.  These  images  are  the  result  of  the  max  iteration  level  set  to  700  but 
they  already  show  improvement  over  figure  30. 
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Coarse  Search  Fine  Search 


(a)  (b) 

Figure  32:  The  algorithm  was  given  enough  iterations  to  converge  at  the  best  estimate 
of  i?o  and  therefore  the  best  estimate  of  the  target  object.  The  max  iteration  level  was 
set  to  1500  but  it  converged  way  before  reaching  this  maximum  number.  Therefore 
the  algorithm  would  also  converge  at  the  same  point  for  this  dataset  at  any  iteration 
level  above  1500. 

4-6  COV  Method’s  Utility  For  Other  Data  Types 

This  research  focused  only  on  simulated  and  laboratory  acquired  electro-optical 
imagery,  but  this  algorithm  should  in  theory  work  on  other  types  of  data.  As  long  as 
the  blurring  function  and  noise  can  be  modeled  it  should  not  matter  if  the  data  is  in 
photons  or  represented  in  the  form  of  spectral  or  radar  frequency  returns  to  name  a 
couple.  As  long  as  the  correct  EM  equations  are  used,  based  on  the  distribution  of 
the  variance  in  the  image,  the  COV  method  should  be  able  to  estimate  the  the  target 
object. 

4-7  Limiting  Factors  That  Effect  Results 

SNR  and  integration  time  and  their  effects  are  analyzed  in  this  section. 

4-7.1  Data  SNR.  The  SNR  of  an  observed  image  is  not  something  that  can 
be  directly  controlled  by  the  image  analyst.  It  is  determined  by  the  the  type  of  target, 
the  energy  source  used  to  measure  it,  and  the  system  which  captures  the  image.  The 
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SNR  of  the  observed  image  drastically  affected  the  results  of  the  COV.  When  the  SNR 
was  too  low  this  method  was  not  always  able  to  determine  the  correct  estimate  of  Rq 
and  it  tended  to  converge  before  the  COV  method  was  able  to  go  through  enough 
iterations  to  remove  the  blurred  effects  of  the  image. 

4-7-2  Integration  Time  Limits.  The  integration  time  or  number  of  iterations, 
on  the  other  hand,  is  up  to  the  image  analyst.  The  COV  method  while  extremely 
accurate  with  a  good  SNR  is  a  time  consuming  process.  The  number  of  iterations 
is  controlled  by  the  user.  As  long  as  time  is  not  an  issue,  allowing  the  algorithm  to 
perform  more  iterations  will  result  in  more  accurate  image  estimations.  Choosing  a 
smaller  number  of  iterations  could  cause  the  algorithm  to  choose  a  larger  Rq  value 
than  it  should  have,  resulting  in  very  badly  blurred  imagery.  Sparse  scenes  can  be 
processed  in  a  matter  of  a  few  minutes  but  full  scene  data  can  take  over  a  hour.  This 
algorithm  in  its  current  form  is  not  appropriate  for  any  type  of  real-time  analysis.  It 
is  better  suited  for  post  analysis  of  imagery. 
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V.  Conclusions  and  Future  Work 


This  section  details  some  conclusions  drawn  from  the  results  of  this  research.  Future 
potential  research  areas  are  also  considered  below. 

5. 1  Conclusions 

First,  has  been  shown  that  the  COV  method  works  on  simulated  data.  It 
was  able  to  effectively  converge  for  the  correct  atmospheric  seeing  parameters  and 
recover  estimates  of  the  true  images,  which  appeared  similar  to  the  actual  images. 
It  was  able  to  perform  effectively  on  fully  illuminated  data  like  both  the  MAP  and 
APEX  algorithms.  More  importantly,  the  COV  method  has  been  shown  to  work  on 
astronomical  data  where  the  other  two  methods  failed. 

Second,  it  was  demonstrated  that  this  method  also  worked  on  laboratory  data. 
This  measured  data  was  generated  in  a  controlled  laboratory  environment.  The  fo¬ 
cusing  device  of  the  imaging  system  was  adjusted  to  produce  a  facsimile  of  our  atmo¬ 
sphere. 

Finally,  the  effects  of  parameters  such  as  SNR,  scene  types,  varying  seeing  pa¬ 
rameters,  and  exposure  time  were  investigated.  These  parameters  were  explored  in  a 
very  cursory  fashion.  There  is  much  work  to  be  done  in  this  area. 

5.2  Future  Work 

There  is  much  exploration  left  in  the  parameters  above.  For  example,  the  effect 
that  SNR  has  on  image  reconstruction  was  studied  but  not  whether  there  exists  a 
certain  threshold  below  which  this  method  will  not  perform  as  expected.  Another  area 
to  look  at  may  be  how  the  COV  method  performs  on  images  that  may  he  somewhere 
between  sparse  and  full-scene.  One  area  which  was  not  researched  was  the  effect  of 
OTFs  with  shapes  that  did  not  £t  the  model. 

It  is  also  important  to  point  out  that  this  technique  was  demonstrated  success¬ 
fully  on  electro-optical  data,  but  by  no  means  is  it  limited  to  this  narrow  band  of  data. 
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It  may  be  possible  to  prove  that  the  COV  method  will  also  converge  on  radar  data, 
microwaves,  x-rays,  etc.  with  the  correct  stopping  criteria  and  EM  update  equations. 


53 


Bibliography 

1.  A.  C.  Likas,  N.  P.  Galatsanos.  “A  Variational  Approach  for  Bayesian  Blind  Image 
Deconvolution”.  IEEE  Transactions  on  Signal  Processing,  Vol.  52,  No.  8:2222- 
2233,  2004. 

2.  A.  MacDonald,  S.  C.  Cain  and  E.  E.  Armstrong.  “Comparison  Of  Registration 
Techniques  For  Speckle  Suppression  In  2-d  Lidar  Image  Sequences” .  SPIE  Confer¬ 
ence  on  Image  Reconstruction  from  Incomplete  Data  III  5558 . 202213.  November 
2004. 

3.  A.  P.  Dempster,  N.  M.  Laird  and  D.  B.  Rubin.  “Maximum  Likelihood  from 
Incomplete  Data  via  the  EM  Algorithm” .  Journal  of  the  Royal  Statistical  Society, 
Vol.  39,  No.  1:1-38,  1977. 

4.  Adam  MacDonald  Lt  Col,  USAF.  Blind  Deconvolution  Of  Anisoplanatic  Im¬ 
ages  Collected  By  A  Partially  Coherent  Imaging  System.  Ph.D.  thesis.  Air  Force 
Institute  Of  Technology,  2006. 

5.  Alfred  S.  Carasso,  David  S.  Bright  and  Andras  E.  Vladar.  “APEX  Method  And 
Real-time  Blind  Deconvolution  Of  Scanning  Electron  Microscope  Imagery” .  Op¬ 
tical  Engineering,  Vol.  41,  No.  10,  October  2002. 

6.  Ce  Liu,  Richard  Szeliski,  William  T.  Freeman  and  Sing  Bing  Kang.  “Noise  Esti¬ 
mation  From  A  Single  Image”.  2006. 

7.  D.  A.  Fish,  A.  M.  Brinicombe  and  E.R.  Pike.  “Blind  Deconvolution  By  Means 
Of  The  Richardson-lucy  Algorithm” .  Journal  of  the  Optical  Society  of  America, 
Vol.  12:5865,  January  1995. 

8.  Goodman,  Joseph  W.  Statistical  Optics.  John  Wiley  &  Sons,  INC.,  2000. 

9.  Leon-Garcia,  Alberto.  Probability,  Statistics,  and  Random  Processes  for  Electrical 
Engineering.  Pearson  Education,  INC.,  2008. 


54 


10.  Loyev,  Vadim  and  Yitzhak  Yitzhaky.  “Initialization  Of  Iterative  Parametric  Al¬ 
gorithms  For  Blind  Deconvolution  Of  Motion-blurred  Images”.  Applied  Optics, 
Vol.  45,  No.  11:2444-2452,  2006. 

11.  R.  Salakhutdinov,  S.  T.  Roweis  and  Z.  Ghahramani.  “Optimization  with  EM  and 
Expectation-Conjugate  Gradient”.  International  Conference  on  Machine  Learn¬ 
ing,  Vol.  20:672-679,  2003. 

12.  X.  Huang,  A.  C.  Madoc  and  A.  D.  Cheetham.  “Image  Multi-Noise  Removal  By 
Wavelet-Based  Bayesian  Estimator”.  Circuits  and  Systems,  Vol.  3:2699  -  2702, 
2005. 


55 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
0MB  No.  0704-0188 


The  public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information,  including 
suggestions  for  reducing  this  burden  to  Department  of  Defense,  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports  (0704—0188),  1215  Jefferson  Davis  Highway, 
Suite  1204,  Arlington,  VA  22202—4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  any  penalty  for  failing  to  comply  with  a  collection 
of  information  if  it  does  not  display  a  currently  valid  0MB  control  number.  PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS. 

1.  REPORT  DATE  (DD-MM-YYYY)\  2.  REPORT  TYPE  3.  DATES  COVERED  (From  —  To) 

24-03-2011  Master’s  Thesis  Sept  2009  —  Mar  2011 

4.  TITLE  AND  SUBTITLE  I  5a.  CONTRACT  NUMBER 


Blind  Deconvolution  Method  Of 
Image  Deblurring  Using 
Convergence  Of  Variance 


5b.  GRANT  NUMBER 


I  5c.  PROGRAM  ELEMENT  NUMBER 


6.  AUTHOR(S) 


I  5d.  PROJECT  NUMBER 


Quentin  MacManus,  Captain,  USAF 


I  5e.  TASK  NUMBER 


I  5f.  WORK  UNIT  NUMBER 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Air  Force  Institute  of  Technology 

Graduate  School  of  Engineering  and  Management  (AFIT/EN) 
2950  Hobson  Way 
WPAFB  OH  45433-7765 

9.  SPONSORING  /  MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

Maj  David  M.  Strong 

Air  Force  Research  Lab,  AFRL/RDSMA 

535  Lipoa  Pkwy. 

Kihei,  Hawaii  96753 

Phone:  808-891-7753 

david .  strong@maui .  afmc .  af .  mil 

12.  DISTRIBUTION  /  AVAILABILITY  STATEMENT 


8.  PERFORMING  ORGANIZATION  REPORT 
NUMBER 


AFIT/GE/ENG/11-26 


10.  SPONSOR/MONITOR’S  ACRONYM(S) 


11.  SPONSOR/MONITOR’S  REPORT 
NUMBER(S) 


APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  UNLIMITED. 

13.  SUPPLEMENTARY  NOTES 

This  material  is  declared  a  work  of  the  U.S.  Government  and  is  not  subject  to  copyright  protection  in  the  United  States. 

14.  ABSTRACT 

Images  are  used  for  both  aerial  and  space  imagery  applications,  including  target  detection  and  tracking.  The  current 
problem  concerning  objects  in  geosynchronous  orbit  is  that  they  are  dim  and  hard  to  resolve  because  of  their  distance. 
This  work  will  further  the  combined  effort  of  AEIT  and  AFRL  to  provide  enhanced  space  situational  awareness  (SSA) 
and  space  surveillance.  SSA  is  critical  in  a  time  when  many  countries  possess  the  technology  to  put  satellites  into  orbit. 
Enhanced  imaging  technology  improves  the  Air  Force’s  ability  to  see  if  foreign  satellites  or  other  space  hardware  are 
operating  in  the  vicinity  of  our  own  assets  at  geosynchronous  orbit.  Image  deblurring  or  denoising  is  a  crucial  part  of 
restoring  images  that  have  been  distorted  either  by  movement  during  the  capture  process,  using  out-of-focus  optics,  or 
atmospheric  turbulence.  The  goal  of  this  work  is  to  develop  a  new  blind  deconvolution  method  for  imaging  objects  at 
geosynchronous  orbit.  It  will  feature  an  expectation  maximization  (EM)  approach  to  iteratively  deblur  an  image  while 
using  the  convergence  of  the  image’s  variance  as  the  stopping  criteria. _ 

15.  SUBJECT  TERMS 


16.  SECURITY  CLASSIFICATION  OF: 
a.  REPORT  I  b.  ABSTRACT!  c.  THIS  PAGE 


17.  LIMITATION  OF 
ABSTRACT 


18.  NUMBER  19a.  NAME  OF  RESPONSIBLE  PERSON 
PAGES  Richard  K.  Martin,  (ENG) 

19b.  TELEPHONE  NUMBER  (include  area  code) 

(937)  255-3636,  x4625;  Richard. martin  afit.edu 


Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std.  Z39.18 


